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P ■ ABSTRACT 

^. 

^ . Context. Modeling of the global heliosphere seeks to investigate the interaction of the solar wind with the partially 

, ^, ' ionized local interstellar medium. Models that treat neutral hydrogen self-consistently and in great detail, together 

with the plasma, but that neglect magnetic fields, constitute a sub-category within global heliospheric models. 

Aims. There are several different modeling strategies used for this sub-category in the literature. Differences and com- 
►^ ■ monalities in the modeling results from different strategies are pointed out. 

w~v ' Methods. Plasma-only models and fully self-consistent models from four research groups, for which the neutral species 

*vj ' is modeled with either one, three, or four fluids, or else kinetically, are run with the same boundary parameters and 

equations. They are compared to each other with respect to the locations of key heliospheric boundary locations and 

with respect to the neutral hydrogen content throughout the heliosphere. 

Results. In many respects, the models' predictions are similar. In particular, the locations of the termination shock agree 
\l ' to within 7% in the nose direction and to within 14% in the downwind direction. The nose locations of the heliopause 

^^ , agree to within 5%. The filtration of neutral hydrogen from the interstellar medium into the inner heliosphere, however, 

OO , is model dependent, as are other neutral results including the hydrogen wall. These differences are closely linked to the 

^^ ' strength of the interstellar bow shock. The comparison also underlines that it is critical to include neutral hydrogen 

into global heliospheric models. 



Key words, heliosphere - interstellar neutral hydrogen - numerical models - kinetic models - multi-fluid models 

1. Introduction to date in the outer heliosphere. Notable sources of in- 

. ... ,. ... formation are the two Voyager spacecraft at a distance 

The mterstellar medmm m the immediate solar neighbor- ^f ^q^ ^^ ^^^ 84 ^U ^^^^ ^^^ g^^ (2007 September), 

hood IS part of the local interstellar cloud (LIC). The respectively, with Voyager 1 having passed i nto the he- 

flow of the partially ionized LIC past the Sun constitutes liogheath region on 2004 December 16 (e.g., iStone et al.l 

a pressure that balances and terminates the expansion of ^^^ ^^^ Voyager 2 on 2007 August 30. For examples of 

the coronal solar wind. These two winds create a morphol- i^-depth analyses of observations relating to the outer heho- 

ogy that includes the termmation shock transition of the ^^^^^.^^ we refer to ot her contribut ions in this special issue 



supersonic solar wind to a hot hehoshcath or hehotail flow. ^ Bzows kilTgl [20081: iPrvor et all [2001 iRichardson et ed] 

An interstellar bow shock is likely to be necessary as well ^2008|; iSlavin fc Frischll2008h . 
to decelerate the LIC flow. The ionized flows of the LIC 

and the solar wind are separated by the heliopause. The Data from the outer heliosphere are sparse, and numeri- 

LIC also supplies the system with interstellar neutrals, pre- cal modeling of the global heliosphere/LIC system plays an 

dominantly with neutral hydrogen (H) . Neutral H interacts important role for the analysis and interpretation of obser- 

weakly with the plasma, mainly through charge exchanges vations. It is needed to relate the undisturbed LIC flow and 

with plasma protons. its physical parameters to the processed and changed flow 

The distance even to the termination shock is large that we observe in the heliosphere inside the termination 

enough that there are only a few in-situ measurements shock. In fundamental ways all the LIC cons traints for- 
mulated in the accompanying papers ([B zowski et al.l 20081: 
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l2008f) involve global heliosphere modeling. Also the eval- 
uation of future data sets from the Interstellar Boundary 
Explorer (IBEX) mission, which focuses on secondary neu- 
trals created in the heliosphere and on the LIC oxygen 
and h elium flow through the heliosphere (jMcComas et alj 
12004) , depends crucially on this kind of modeling. 

All such global models make assumptions and simplifi- 
cations, most often with the goal of isolating the influence 
of a specific physical effect (e.g., the tilt of the LIC mag- 
netic field with respect to the LIC flow vector) , or in order 
to keep computation times reasonable. The identification 
of heliosphcric a symmetries with respect to the helium LIC 
flow vector (Mobius et al.ll2004l : iLallement et al.ll2005l . and 
references therein) has increased interest in the develop- 
ment of realistic, three-dimensional (3D) MHD models, as 
different orientations and strengths of the interstellar mag- 
netic field can help to explain these asymmetries. However, 
the fact remains that neutral interstellar H entering the he- 
liosphere has a more decisive influence on the heliospheric 
shape, extent, and particle content. For this reason, we fo- 
cus here on numerical models that treat the plasma/neutral 
interaction in a self-consistent way, but neglect the influence 
of interplanetary or interstellar magnetic fields. The models 
are, in principle, 3D plasma/neutral codes for which plasma 
and neutrals are coupled by charge exchange. Wherever it 
is possible, the assumption of azimuthal symmetry reduces 
the numerical methods effectively to 2D while still calcu- 
lating the 3D heliosphere. The results of our investigation 
will be also applicable to 3D MHD models (for a recent 
overview, see iPogorelov et al] ()2008f )l. as long as the lat- 
ter also include neutrals self-consistently, as is essential for 
models of the global heliosphere. 

The charge exchange interaction is weak enough that 
the mean free path lengths of neutral H are often large com- 
pared to typical heliospheric distances (see the discussion 
in section 3.2). Neutral H is thought to be in local thermo- 
dynamic equilibrium in the LIC, but as charge exchange 
proceeds in the heliosphere, secondary neutrals arise, and 
they also can exchange charge with plasma ions. This effec- 
tively drives neutrals out of equilibrium in the heliosphere, 
and plasma and neutrals equilibrate again only far away 
from the heliosphere. 

In H-p charge exchange, the newly born (secondary) 
neutral H has the velocity characteristics of the plasma 
protons at the location of interchange. However, the neu- 
tral is no longer bound to the plasma flow and follows a 
simpler trajectory than the underlying plasma parcel. Due 
to this, it is convenient to sort the neutrals into different 
populations depending on their origin. We will label here 
as component 1 the primary neutrals directly from the ISM 
as well as those born in charge exchange outside the bow 
shock. Component 2 are those secondary neutrals that are 
created by charge exchange between bow shock and he- 
liopause. They reflect the conditions of the warmer inter- 
stellar plasma decelerated in the bow shock. Because of the 
deceleration, there is a neutral density increase downstream 
of the bow shock, the hydrogen wall. 

We label as component 3 neutrals those that are born 
from the hot heliosheath and heliotail plasma, and compo- 
nent 4 those born in the supersonic solar wind between the 
Sun and the termination shock. Component 3 neutral ve- 
locities are dominated by the large thermal proton velocity 
of the heliosheath and heliotail, and hence their direction 
is mostly random. Since a fraction of component 3 neutrals 



are directed to the innermost heliosphere and can be de- 
tected as energetic neutral H, this whole component 3 is 
often referred to as "heliospheric ENA" (energetic neutral 
atoms). The fourth neutral component has recently been 
called "neutral solar wind" (NSW) because its cold, fast 
velocity characteristics are similar to those of the super- 
sonic wind of the inner heliosphere. Note that in spite of its 
name, the NSW as defined here is distinct from the neu- 
tral hydrogen originati ng from the Sun fe.g.. lBlum fc Fahii 
llQTOtlOlsen et al.lll994^ ■ which has been called neutral solar 
wind earlier as well. 

In the non-MHD models that are currently applied to 
the global heliosphere problem there is agreement that due 
to the out-of-equilibrium nature of neutral H it needs to be 
modeled separately from the ionized matter. The plasma 
is commonly modeled by gas-dynamic methods. There are 
two different popular methods for treating the neutrals, to 
be coupled to the plasma in a self-consistent way. The first 
method is kinetic, where particle methods such as Monte- 
Carlo simulate the neutral populations on a Boltzmann- 
microscopic level. The kinetic treatment is motivated by 
the usually large mean free paths of neutrals. The sec- 
ond approach, the multi-fluid method, is to simulate each 
of the four neutral components as a separate fluid on an 
Euler-macroscopic level, and assumes that the superposi- 
tion of the resulting four Maxwellian distributions repre- 
sents the true, generalized distribution function of helio- 
spheric neutral H well. Sometimes, fluid models are being 
restricted further (by choice) by decreasing the num ber of 
fluids to less than four, as in the IZank et all (|l996l ) multi- 
fluid mo del (compo n ent I and 2 combined into one fluid) 
and the iFahr et aLl (120001 ) Bonn model (components 1-4 
combined, but fluids describing pickup ions and cosmic rays 
introduced) . 

Without going into any detail, the two neutral model- 
ing approaches can be summarized as follows. The main 
advantage of the kinetic approach is that it does not re- 
strict the shape of the neutral distribution function, and 
thus allows the irregularity of the neutral distribution in 
the heliosphere to persist everywhere in the heliosphere. 
The main advantage of the fluid approaches is that they 
are orders of magnitude faster computationally, and that 
their usual field variables (density, velocity, pressure) are 
smooth down to the grid and timestep resolution. The main 
disadvantage of particle kinetic methods is that their ac- 
curacy is driven by particle statistics, i.e. to increase the 
accuracy of results for a particular location at a particu- 
lar time, more particles have to be generated to coincide 
there at the desired time. Both kinetic models below (sec- 
tion 3) employ variable particle weights to allow trajec- 
tories to be split, leading to a significant improvement in 
the statistical accuracy (to a targeted ~2% level) at rea- 
sonable com£utationaJ_costs. The splitting procedure used 
in the Baranqv fc Malarnal (11993) model is descri bed by 
iMalamal (|l99lh . while the lHeerikhuisen et all (|2006[ ) model 
uses a similar method based on splitting during charge ex- 
change. The main disadvantage of fluid models is that each 
neutral component is forced into local thermal equilibrium, 
which, at the very least, constitutes a loss of information 
(Maxwellians instead of a more general distribution). 

Two studies have recently engaged in detailed compar- 
isons of kinetic models versus multi-fluid gasdynamic mod- 
els of neutral hydrogen in the heliosphere, and put for- 
ward some of the possible physical reasons for the differ- 
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ences that invariably occ ur (jAlexashov fc Izm.odenovll2005l : 
iHeerikhuisen et al] |2006') . For their study comparing global 
hehospheric models, (Alcxashov & Izmodenov (2005, here- 
inafter AI05) set out using the Moscow kinetic code devel- 
oped over the years in Moscow starting from the original 
iBaranov et alj (|l99lh kinetic-gasdynamic model. A certain 
set of solar wind and interstellar boundary parameters is 
used throughout their study. AI05 compare the kinetic re- 
sult, and non-self consistent variants of it, to multi-fluid 
models in which neutral H is modeled by one to four fluids, 
coupled self-consistently to the same gas-dynamic plasma 
code used also for the kinetic model. They find that ki- 
netic and multi-fluid models never agree completely. The 
agreement with the kinetic method is best for the four-fluid 
model, and worst for the one-fluid model. The boundary 
locations are further out in the four-fluid case when com- 
pared to the kinetic model, namely, the upwind bow shock 
(BS) by 4%, the termination shock (TS) by 5%, and the he- 
liopause (HP) by 9%. The hydrogen wall material is more 
decelerated and consequently the peak density is larger, 
and the filtration more severe (less H entering through the 
termination shock). This discrepancy is started by the ki- 
netic model having a weaker (plasma) bow shock than the 
four-fluid model, which in turn is likely caused by more sec- 
ondary neutrals passing to the interstellar side of the bow 
shock than in the fluid case, as displayed by AI05. 

In a similar investigation, IHeerikhuisen et all (|2006l 
hereinafter HFZ06) use their own kinetic code and com- 
pare the result with their own version of a four-fluid model. 
They, too, find a weaker bow shock in the kinetic case 
compared with the fluid case, and the same chain of con- 
sequences, including a larger H density in the hydrogen 
wall and a smaller density passing through the termination 
shock for the fluid case. While TS and HP are farther out 
for the fluid case, the BS is less far than in their kinetic 
model. At four representative locations in the heliosphere, 
HFZ06 compare the parallel velocity distribution function 
between kinetic and multi-fluid models, the latter being 
a superposition of 4 individual Maxwellians. At least for 
those 4 points on the stagnation axis, they flnd that the 
two interstellar neutral components coincide very well with 
the Maxwellians from a four-fluid model, and only the two 
heliospheric neutral components deviate. AI05 and HFZ06 
agree that the NSW component is much hotter (i.e., broader 
velocity distributions) in the kinetic model than in the fluid 
model. The EN A- "fluid" is the most problematic to be fit 
by a Maxwellian, at least outside the inner heliosheath and 
the hcliotail. 

HFZ06 also compare their results to those obtained by 
AI05 with their codes. This comparison is possible because 
HFZ06 use the same boundary parameters. The two four- 
fluid codes correspond very well to each other, minus a 
subtle difference that comes about by the different internal 
treatment of the bow shock in the two plasma codes. The 
two kinetic code results also differ in hydrogen density be- 
tween bow shock and heliopause, which again might have 
to do with the internal treatment of the bow shock in the 
underlying plasma codes. 

In this paper, we take one step further back and com- 
pare the results of sophisticated, comparable global helio- 
sphere models (albeit all axisymmetric and non-MHD), run 
on the same boundary data set characterizing the solar 
wind at 1 AU and the pristine interstellar medium. Since 
the modeling strategies even for the plasma gas-dynamic 



model are different across the four groups considered (they 
are compared in section 2), the differences between the 
models are going to be larger than for the case of the in- 
ternal comparisons of AI05 and HFZ06. In this sense, the 
paper focuses not on discovering additional physical rea- 
sons for differences between the kinetic and the multi-fluid 
approach. Rather, we are trying to state quantitatively how 
far apart or close some key results are, in order to give the 
wider community a sense how accurate statements derived 
from current neutral/plasma models likely are. 

The models used within the observational contributions 
of this special issue all are related to the global models 
outlined below, and all make use of specific additions or 
modifications depending on the specific issues addressed. 
The Richardson and W ang one-dimensional MHD model 
(jRichardson et al]l2008^ concentrates on the solar wind - 
interstellar flow interaction to describe the slowdown in the 
supersonic solar wind. It does not include many of the intri- 
cacies of the global models beyond the termination shock as 
discussed below, but incorporates the detailed solar wind 
temporal structure to come to a m eaningful comparison 



iparison 
T (|2n08t ) 



between inner and outer heliosphere. iBzowski et al 
take the interstellar flow from a kinetic model similar to 
the one in AI05, and then add a Monte Carlo calculation 
of the history of individual particle trajectories in the inner 
heliosphere to g et the picku p ion c haracteristics between 1 
and 5 AU. For IPrvor et all (|2008[ ) it is important to add 



the radiation transport of solar Ly-a, including multiple 
scattering on a global heliospheric model. In this sense, the 
global models discussed below may serve as proxies for the 
modeling used for the entire special section. 

The remainder of the paper is organized as follows. In 
section 2, we compare results from single fluid, plasma-only 
models, and in section 3, from flve neutral/plasma global 
heliospheric models for one particular parameter set. In sec- 
tion 4, we attempt to qualify our results and put them into 
perspective of the overall goal of realistic models of the 
global heliosphere. 

2. Comparison of plasma-only models 

The five self-consistent global heliospheric models that are 
compared in this paper are the lBaranov fc Malama' (1993) 
M oscow model ("BM " ), th e IGPP-UCR kinetic model 
bv IHeerikhuisen et all (I2006D ("H ee") and the multi-fluid 
mo del bv iFlorins ki et al.' ('2005'^ ("Flo") ext ended from 
the iFlorinski et a l. (2003) two-fluid model, th e iPauls et al. 



(TOOS*) style multi-fluid model modified by 'Miille r et al 



(20061 ("Mue"), and the Bonn five-fluid model (Fahr et al. 



I2000h as used bv lScherer fc Fahil (|2003^ ("Sch"). Ah these 
models use a gas-dynamic description of the plasma, and 
therefore we start with a comparison of the plasma-only 
part, i.e. the Euler equation solvers used by these groups for 
their plasma part under the assumption that there are no 
source terms on the right-hand-sides of the fluid equations 
(no neutrals in the system). All groups ran their plasma 
code with the solar wind and interstellar boundary con- 
ditions listed in Table 1. It was assumed that the plasma 
consists of equal (comoving) densities of protons and elec- 
trons, in other words, the thermal plasma pressure equals 
twice the thermal proton pressure. The magnetic field, as 
well as the solar gravity, were neglected. 

As can be expected, the results from the four groups are 
very close to each other (the plasma parts of the Hee and 
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Fig. 1. (a) Number density profiles for the plasma-only models of all four groups, in the directions upwind (0), crosswind 
(90), and downwind (180) with respect to the LIC flow, respectively, (b) Detail around upwind BS. 



Table 1. Boundary parameters, plasma-only models. 



Variable 



lAU LIC [units] 



proton density 

velocity 

temperature 



7 0.06 [cm-^] 
375 26.4 [kms-^] 
73,640 6530 [K] 



Table 2. Key results from plasma-only models. 



Result 


BM 


Flo 


Mue 


Sch 


mean 


upwind TS [AU] 
upwind HP [AU] 
upwind BS [AU] 
downwind TS 


149 
196 
340 
385 


148 
199 
340 
319 


153 
204 
365 
396 


152 
212 
380 
304 


150.5 ± 2.4 

202.8 ± 7.0 

356 ± 18 

351 ± 46 



the Flo models are identical, and not hsted separately). 
Figure [T] shows this with the number density profiles of 
all four models, in three representative directions in the 
commonly used heliocentric reference frame. In all three 
directions, all the densities follow a r~^ power law in the 
supersonic solar wind before encountering the termination 
shock. The upwind termination shock (TS) distances are 
very close to each other, whereas downwind there is more 
variability. The shock strengths (ratio of downstream to 
upstream density) are more or less the same. Also the den- 
sity contrast across the heliopause (HP) is similar across 
the four models. Also obvious is that the BM model uses 
capturing methods to identify and enforce discontinuities, 
whereas the other three models do not employ such tech- 
niques, and transitions are spread over a few grid points 
(see the bow shock of the Mue model for an example, Fig. 

Table [2] lists some key results for the shock and he- 
liopause locations. The similarities in the results are evi- 
dent. The last column comprises a simple average across 
the four models for each result, and the standard deviation 
hints at the range that the results span. The different mod- 
els basically agree on the upwind TS and HP locations, and 
are a little bit more spread for the BS, and yet more for the 
downwind TS results. 

Besides the treatment of shocks and discontinuities, 
there are obviously many other reasons why the four mod- 
els vary from each other. Each of the four models makes 
different choices related to the grid configuration, resolu- 



tion, and the extent of the simulation domain. Also, there 
are four different choices of the numerical transport and 
diffusion schemes to solve the Euler equations. BM use a 
Godunov-type numerical scheme with moving adaptive grid 
while capturing three discontinuities — the heliopause as 
contact discontinuity, and the termination and bow shocks. 
The accuracy of the numerical scheme resolution is im- 
proved by using a "minmod" limiter. The plasma part 
of the numerical algorithm of the Flo multi-fluid model 
uses the total variation diminishing (TVD) finite volume 
scheme based on the Courant-Isaacson-Rees approximate 
Riemann solver. Conservation laws for the neutral compo- 
nents (section 3) are solved using the more diffusive TVD 
Lax-Friedrichs method. The Hee plasma part is identical to 
that of Flo. The ZEUS-3D algorithm underlying the Mue 
model is based on the method of finite differences on a 
staggered mesh, incorporating a van Leer monotonic ad- 
vection scheme, and von Neumann-Richtmyer artificial vis- 
cosity at shock fronts. For all fluids in the Sch model the 
Euler equations are formulated for quantities conserving 
the flux of mass, momentum, and energy, and are subjected 
to second order Riemann solvers using the Lax-Friedrichs 
method with an entropy fix. For large pressure gradients a 
Harten-Lax-van Leer solver is implemented. 

The high-Mach number regime of the supersonic solar 
wind is an instructive example of the modeling technique 
differences, and their consequence for heliospheric studies. 
While each technique is optimized to conserve crucial quan- 
tities (for example, mass flux from grid cell to grid cell), the 
calculations of density, velocity and pressure deviate from 
model to model. Small flux errors are evident in Figure [5] 
showing the conserved total particle flux nvr'^, where the 
ideal value (2.625 x 10^ AU^ cm^^ s~^) is approximated well 
by the BM model. The Flo model also conserves this quan- 
tity, albeit a smaller value was introduced at the bound- 
ary. Immediately upstream of the TS, the modeled densi- 
ties range from 6.75 to 7.35 cm~^/r^, and the velocities 
(ideally 375 km s~^) range from 376 to 383 km s^^. As the 
location of the TS is determined by the ram pressure at 
the TS balancing the ISM pressure, the subtle variations in 
ram pressure are a natural explanation of the TS differences 
in Table [21 Similar effects explain the other discrepancies. 
The BS distances basically follow the trend of the HP dis- 
tances, as the BS shock compression ratios are quite similar 
between all four models (2.2-2.3; cf. Fig. [Us). We note in 
passing that the stagnation axis from which data for Table 
[2] and Figures [1] (except for 90°) and[2]are taken, is numer- 
ically somewhat problematic in that for the axisymmetric 
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Table 3. Boundary parameters, full models. 
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Fig. 2. Upwind profiles of flux conservation nvr"^ in the 
supersonic solar wind. 



models at hand, it actually consists of boundary grid zones 
and not of interior zones. 

One specific difficulty faced by every model of the he- 
liosphere is the treatment of outflow boundary conditions 
in the heliotail. Because the tail plasma flow is subsonic, 
the boundary is influenced by waves and disturbances en- 
tering from downstream, i.e., from the regions not included 
in the simulation. Simple outflow boundary conditions are 
used by the Mue and Sch models, in effect copying inte- 
rior solutions to boundary shadow zones and thus making 
derivatives at the boundary small. This approach could lead 
to waves reflected off the boundary and reentering the sim- 
ulation, an unphysical situation. A somewhat more compli- 
cated approach, employed by the Flo model, is to apply a 
"non-reflective" outflow boundary condition, whereby the 
flow is reaccelerated to a sonic point through an insertion 
of a rarefaction fan at the boundary of the domain. In the 
BM model the tail computation region is extended up to 
the region where the solar wind is supersonic again. While 
the other models agreed to an outer boundary at 1000 
AU, the BM model extends the simulation domain to 6000 
AU tailward for this rea son (jlzmodenov k, Alexashovll2003t 
lAleksashov et al.l 12004 ) . at the expense of resolution and 
computation time. Regardless of their degree of sophistica- 
tion, it should be realized that all tail boundary conditions 
used by heliospheric models are not physically exact in the 
strict sense, except for the BM model where the outer tail 
flow is supersonic and, therefore, the boundary conditions 
are correct. The boundary handling contributes to notice- 
able differences in the distances of the TS in the downwind 
direction predicted by different models (see Table [2]), yet it 
is not the only issue involved, as there is a curious pairing 
of BM and Mue model distances on the the one hand, and 
Flo and Sch models on the other hand. 



3. Comparison of self-consistent plasma/neutral 
models 

3.1. Model results 

We now proceed to introduce neutral interstellar hydrogen 
(H) into the system and, using the plasma codes of sec- 
tion 2, switch on the full, self-consistent plasma/neutral 
codes in which the plasma and neutral H influence each 
other through appropriate source terms. All groups calcu- 
late their global heliosphere with the solar wind and inter- 
stellar boundary conditions listed in Table [31 Again, the 
magnetic fields are neglected, as are gravity and radiation 
pressure. The H-p charge exchange cross section depends 



Variable 



lAU Lie [units] 



proton density 
H density 
velocity 
temperature 



7 0.06 [cm-3] 

- 0.18 [cm-3] 

375 26.4 [kms-i] 

73,640 6530 [K] 



on the relative velocit y; for this paper all five models use 
the iMaher fc Tinsley (jl977t ) cross section. A photoioniza- 
tion rate of 10~* s~^ (1 AU / r)^ is assumed, and other 
ionization channels such as electron impact ionization are 
neglected throughout. 

Figure [3] gives a good overview of the results for the two 
available plasma - kinetic neutral models BM and Hee, and 
the three multi-fluid models Flo, Mue, and Sch. The left 
panel shows the plasma density profiles for upwind (solid) , 
crosswind (dotted), and downwind (dashed) directions, and 
the right panel displays the information for total neutral 
density in the same format. The plasma results exhibit 
a level of similarity to each other that is comparable or 
slightly better than the level of similarity of the plasma 
results in the previous section (Fig. [T]). In particular, the 
upwind HP location nearly coincides in all five models. The 
first four entries of Table [3] contain the key locations of the 
heliosphere for the full models. Again, the standard devi- 
ations in the last column express the range of the results 
against a simple arithmetic mean of the values in each row. 
It can be seen that also the upwind TS locations agree to 
within 7%, and only BS and downwind TS disagree (up to 
14% each). While this type of disagreement was also found 
in the plasma-only cases of section 2, it tends to now be 
larger, especially for the bow shock. We note in passing the 
dramatic effect on the heliosphere boundary locations that 
the inclusion of neutrals has. The results of Table [5] are 
significantly larger than those of Table 2] even though the 
boundary parameters of the plasma-only case are identical 
to those of the plasma/neutral case. 

In the neutral H density (Fig. [31 right) all models ex- 
hibit an overdensity (hydrogen wall) downstream of the bow 
shock, and a subsequent rapid drop in the density approach- 
ing the heliopause and further inside. For this and simi- 
lar diagnostics, the neutral multi-fluid results are summed 
(averaged) into a single total neutral hydrogen quantity, 
in the simplest manner as ritot = Xli "-« ^^^ total density, 
vtot = n:irJi. ■ J2i ^iVi for velocity, and T = n^j • Y,i ntTi for 
temperature. 

The two fluid models with less than four neutral flu- 
ids (Mue, Sch) almost agree in the sharpness and the 
peak height of t he hy drogen wall. As in previous findings 
(jBaranov et all (|1998[ ): Figure 2 bv iMcNuttI (|2004D : AI05; 
HFZ06), the hydrogen wall is quite a bit higher for these two 
fluid models compared to the kinetic models BM and Hee. 
The two latter models match each other well in neutral hy- 
drogen. The hydrogen wall of five-fiuid model Flo is higher 
than the kinetic ones, but closer to those than to the other 
multi-fluid models. The peak densities in the hydrogen wall 
are listed in Tabled The hydrogen wall profiles fit the gen- 
eral trend displayed in Figure 3 of AI05: There, the one- 
fiuid model (most similar to the Sch model) resulted in the 
highest-peaked hydrogen wall, the three-fiuid model (most 
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Fig. 3. Density profiles of plasma (left), and of neutral H (right), for three directions, for the five self-consistent 
plasma/neutral models. 

Table 4. Key results from plasma/neutral models. 



Result 


BM 


Hee 


Flo 


Mue 


Sch 


mean 


upwind TS [AU] 


87 


85 


90 


94 


96 


90.4 ± 4.6 


upwind HP [AU] 


130 


126 


132 


130 


138 


131.2 ±4.4 


upwind BS [AU] 


245 


274 


260 


236 


209 


245 ± 25 


downwind TS 


177 


166 


190 


214 


192 


188 ± 18 


BS compression ratio 


1.2 


1.2 


1.4 


1.7 


2.3 


1.6 ±0.5 


peak uh [cm^ ] 


0.27 


0.27 


0.30 


0.43 


0.46 


0.346 ± 0.092 


UH at TS [cm-^] 


0.134 


0.125 


0.109 


0.094 


0.126 


0.118 ±0.016 


filtration / 


0.74 


0.69 


0.61 


0.52 


0.70 


0.65 ±0.08 


vh at TS [km s"^] 


20.7 


20.8 


23.4 


21.3 


19.2 


21.1 ±1.5 


Th at TS [K] 


26800 


30900 


15500 


21000 


12000 


21200 ± 7800 



similar to the Mue model) exhibited a somewhat smaller hy- 
drogen wall with a very sharp rise on the interstellar side, 
and the BM model had a small peak, with a smooth H den- 
sity rise and fall, that was relatively closely matched by a 
four-fluid model (most similar to the Flo model). Note that 
the neutral H column density through the upwind direction 
is basically constant; the displayed different hydrogen walls 
are either tall and narrow, or smaller and broad. 

The largest contributor to the differences between the 
simulated hydrogen walls is the distribution of plasma ve- 
locities, notably the component parallel to the ISM flow. 
FigureSlshows an overview of the radial velocity component 
as a function of heliocentric distance in the left panel, and 
the right panel zooms in around the bow shock distance. 
The two kinetic models display a plasma deceleration up- 
stream of the bow shock due to charge exchange with com- 
ponent 3 and 4 neutrals that are streaming antisunward and 
have passed the bow shock. To a lesser extent, the Flo and 
Mue models exhibit a similar deceleration, whereas the Sch 
model cannot resolve such counterstreaming fluid elements 
(they deposit their energy already far downstream of the 
BS). As a consequence of the deceleration both upstream 
and downstream of the BS, the BS is the weakest in the 



kinetic cases, followed in shock strength by the Flo model, 
and is the strongest in the one-fluid case, with the Mue 
model in between (Table |31). The shock-capturing method 
of BM arrives at a very weak BS. This range of bow shock 
strengths explains the more gradual hydrogen wall in the 
kinetic cases. The hydrogen wall is of lesser amplitude in 
the kinetic cases because the velocities downstream of the 
BS are distinctly larger (absolute magnitude) than those 
in the fluid cases, and therefore charge- exchanged neutrals 
are not decelerated as much as in the one-fluid case where 
the plasma velocity is the smallest. To appreciate this dif- 
ference, one has to mentally shift the plots of Fig. 0] (right) 
so that the individual bow shocks line up. As expected, a 
stronger bow shock results in a more decelerated plasma, 
and therefore a larger peak density of the hydrogen wall. 
Typically, this also means a lesser distance of the BS to the 
Sun, and this trend is evident in Tabled] The exception is 
model BM that experiences additional deceleration down- 
stream of the BS, such that the BS standoff distance is not 
as far outward as the shock strength would suggest. 

The different hydrogen walls result in different neutral 
densities downwind of the TS, where the neutrals enter re- 
gion 4 of the supersonic solar wind. The filtration ratio 
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Fig. 4. Plasma radial velocity profiles of all models for three directions (left), and a magnification around the upwind 
bow shock (right). 



nH{TS)/nH{oo) is listed in Tabled along with the two 
absolute densities discussed above. In principle, it can be 
expected that the slower the outer heliosphere plasma is 
(i.e. the higher the hydrogen wall), the more neutrals get 
deflected around region 4 and hence the stronger the fil- 
tration is. However, for the slowest case, the one-fluid Sch 
model, the small width of the inner heliosheath compen- 
sates for this effect and cr eates TS neutral densitie s similar 
to the kinetic BM model. iRichardson et all (|2008f ) investi- 
gate the slowdown of the supersonic solar wind due to the 
pickup process derived from neutrals. They calculate a 15% 
solar wind slowdown between 5 AU and 78 AU from a 2006 
conjunction of the Ulysses and Voyager 2 spacecraft. The 
five models described here indicate less of a slowdown in 
this distance range, namely, between 7% (Mue, Sch) and 
10% (aU others). 

While the differences in the neutral density at the ter- 
mination shock (and consequently, the differences in the 
filtration ratio) are in the 20% range, the predicted neutral 
velocity at the termination shock varies much less among 
the five models. Table |4] contains the corresponding results, 
as well as those referring to the (total) neutral temperature 
at the termination shock. The velocities are within 6% of 
each other; in other words, all five models predict a slow- 
down of the neutrals by (5.3 ± 1.5) km s~^ because of their 
passage from the LIC to region 4. The picture is much less 
clear for the neutral temperature, where large variations 
exist between the two-fluid model (smallest T) on the one 
side, and the kinetic models (largest kinetic T) on the other 
side, with the Mue model in between. 

3.2. Discussion 

As is evident from the results above, the two kinetic mod- 
els yield very similar results, while the discrepancies be- 
tween the kinetic models and the one-fluid model are the 
largest, which was also found by AI05. These two model- 
ing alternatives bracket the range within which the rest of 
the models fall. The kinetic codes have no provision for di- 



rect neutral-neutral interaction, and hence have no direct, 
built-in mechanism that would drive the neutral distribu- 
tions toward Maxwellian equilibrium. On the other hand, 
the Euler equations of gas dynamics are derived under the 
assumption of frequent particle collisions (even though they 
seem to be valid beyond that range) to yield thermody- 
namic equilibrium, and hence the single-fluid description 
of neutral H is related to envisioning frequent neutral col- 
lisions. The multi-fluid codes are interesting in that they 
disallow neutral-neutral interactions between neutrals that 
are created in different thermodynamic regimes (which usu- 
ally means that they are distinctly separated in velocity 
space) , while still envisioning neutral collisions within those 
regimes themselves. As neutral H is being modeled with 
more and more fluids, fewer and fewer neutral collisions 
are assumed, which results in a convergence toward the ki- 
netic results, seen in the results above and those by AI05 
and HFZ06. 

The absence of neutral-neutral interactions, in par- 
ticular in the kinetic models, does however not neces- 
sarily mean that the distributions are completely non- 
Maxwellian. Charge exchange with plasma protons injects 
neutrals drawn from a Maxwellian distribution, hence the 
secondary neutrals have a tendency to organize in distri- 
butions similar to Maxwellians. In this sense, charge ex- 
change constitutes an indirect (and inefficient) channel for 
neutral equilibration. The charge exchange mean free path 
(nifp) is nowhere very small in the heliosphere, but is some- 
times small enough for charge exchange to occur frequently, 
driving neutrals toward equilibrium. Indeed, some example 
mfps derived from the Mue model are quite short. In the 
region upstream of the bow shock, the interstellar neutral 
mfp is ~ 200 AU, and the mfp of neutrals having been 
generated in regions 3 and 4 (the regions occupied by solar 
wind) and having streamed to upstream of the bow shock is 
even smaller, ~ 100 AU. For the outer heliosheath (between 
BS and HP), mfps are below 100 AU, and go down even to 
20 AU close to the nose of the heliopause (on the interstel- 
lar side) where plasma velocities become small. Therefore, 
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the Knudsen number in these regions is small, down to 
~ 0.2. In these instances, and even when Knudsen numbers 
are larger, forcing the neutrals into multi-component fluids 
with Maxwellian distributions in general does not change 
their distribution much. This interpretation is backed up by 
the findings by HFZ06 who decompose the neutral distribu- 
tion function into the contributions from the four compo- 
nents. They find that the interstellar component (compo- 
nent 1) and the outer sheath component 2 do behave like 
Maxwellians even when treated fully kinetically (cf. their 
figure 5). In this context, AI05 argue that one of the more 
fundamental differences between kinetic and fluid treat- 
ments of neutral H is related to the fact that for the in- 
terstellar temperature of 6530 K, the thermal velocity of H 
is already about half of the bulk velocity value of 26 km/s. 
This means that individual ISM particle trajectories have 
sizeable perpendicular velocity components that are repre- 
sented in the kinetically modeled trajectories. In contrast, 
the fluid description of the same region has a strictly par- 
allel bulk velocity, and effects of the perpendicular particle 
motion are handled by a non-zero neutral pressure. 

For secondary neutrals produced in the solar wind (re- 
gions 3 and 4), the fluid picture is capturing some aspects 
of the particle behavior less well. Component 3 might be 
reasonably approximated as a hot Maxwellian in region 3, 
but component 3 streaming out of this region will have a 
complicated distribution function. For points outside the 
HP, the distribution function will be resembling a half- 
Maxwellian, with velocity components toward the HP miss- 
ing (half-Maxwellian in region 1, HFZ06). For region 4, 
where component 3 neutrals constitute the energetic neu- 
tral atom (ENA) hydrogen background, the distribution 
is complex (each location is reached in principle by ENA 
from all heliosheath and heliotail positions), and certainly 
non-Maxwellian. The fluid approach consequently has com- 
ponent 3 very hot and with a small velocity in region 4. 
Similar findings apply to component 4 neutrals, which are 
cold and fast in principle. As they stream to distant loca- 
tions in region 1, the kinetic codes allow their distribution 
function to broaden unhindered by interactions and thereby 
gaining a large kinetic temperature, whereas the fiuid com- 
ponent 4 experiences the adiabatic cooling of the radial 
expansion, and ends up with much smaller temperatures. 
The differences in component 3 and 4 neutrals present in 
region 1 between the kinetic and the fiuid picture set the 
stage for the different bow shock strengths mentioned above 
and hence infiuence the BS location and the hydrogen wall. 
In the solar wind region, the absolute energy transfer to 
the plasma due to charge exchange by component 3 and 
4 neutrals seems insensitive to the subtle differences in the 
distribution function there, and hence TS and HP locations 
are basically unaffected. 

4. Sources of Error 

Using a multi-fluid approach instead of a particle kinetic 
method incurs a systematic error in the neutral distribu- 
tions, and therefore also in the plasma distributions. This 
has been discussed in the previous section and in the lit- 
erature (e.g., AI05 and HFZ06). In this section we want 
to discuss the source of other systematic errors that con- 
tribute as well to differences between any modeled global 
heliosphere and the real system as observed through helio- 
spheric measurements. 



4.1. Numerics 

As illustrated in section 2, simple choices for the funda- 
mental algorithm for following the non-MHD fluid equa- 
tions, combined with choices for grid resolution and orga- 
nization (e.g., spherical vs. Cartesian, or fixed resolution 
vs. location-dependent) and choices relating to the extent 
of the computation domain determine the outcome of even 
the simplest, plasma-only heliospheric simulation. Different 
choices will conserve different quantities better, usually at 
the expense of other quantities (see the conservation of nvr'^ 
in Fig. ID). 

Another common issue of fiuid simulations is the han- 
dling of discontinuities such as termination shock, bow 
shock, and hcliopause. The solutions used by the models 
in section 2 range from smearing out discontinuities over 
three grid cells to shock capturing methods that supply the 
discontinuous solution externally, and not from the fluid 
algorithm used everywhere else. The treatment of disconti- 
nuities is hence quite sensitive to the local grid resolution 
at the heliospheric boundaries. 

As is usually the case, both multi-fluid simulations as 
well as billion-particle simulations involve a myriad of in- 
dividual steps, with the potential that even numerical ac- 
curacy comes into play as a potential source of error. This 
is presumably less important, however, as the simulations 
eventually settle into a converged, time-independent state 
for which roundoff errors should cancel. 



4.2. Cross sections 

The results of global heliosphere modeling are sensitive to 
the cross sections that are chosen for the resonant charge 
exchange between protons and neutral hydrogen. Often, 
studies of this charge exchange cross section have been mo- 
tivated by investigations of the terrestrial ionosphere inter- 
acting with the terrestrial neutral exosphere, and hence are 
not meant for higher energies > 1 keV, which is a source 
of error for charge exchange involving component 3 and 4 
ne utrals. 

llzmodenov et al.l ()2000f ) and iFahr et all (|2007[ ) have 
reviewed issues relating to the relevant cross sections. 
Typically heliospheric modelers have adhered to the energy- 
dependent cross sections by Maher & Tinsley (1977) and 
iFite et ahl (jl962l ). Both are fitting formulae of the form 
(a — blogv)^, where v is the relative velocity be- 
tween the interaction part ners. A recent compilation by 
iLindsev fc Stebbing^ (|2005[ ) arrives at a yet different cross 
section approximation. The cross sections are still uncertain 
to approximately 10% (solar wind speeds) and up to 40% 
(slow speeds; see, e.g.. lBzowski et al] (|2008[ )). not only be- 
cause of the fitting itself and the extrapolation of these fits 
beyond their intended velocity range, but also the underly- 
ing experimental data from different groups do not always 
reconcile easily. 

Heliospheric modeling is very sensitive to the actual 
cross section values. In order to not repeat work reported 
elsewhere, we would li ke to draw attention to figure 4 by 
HFZ06 and figure 8 bv lBaranov et all ()1998[ ). Each of these 
two figures compares two (respective) heliospheric models 
that differ only by the choic e of the cross section, i.e. ei- 
th er using th e value s by Mah er fc TinslevI (1977) or those 
bv lFite et al.l (|l962f ). For both papers, the results indicate 
that the shift in heliospheric boundary locations like TS 
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and HP is minor, about 1-3%, but that the consequences 
for the hydrogen wall and for the neutral filtration factor are 
quite a bit larger : The hydrogen wall is ~ 14% larger in the 
iFite et all (|1962D case, and correspondingly, there is more 
filtration going on for that case (sm aller / number, a s used 
in Table [J). The reason is that the iFite et all (I1962D cross 
section is larger than the lMaher fc TinslevI ( 1977f ) cross sec- 
tion at key energies. 

In this sense, the uncertainty in the cross section and 
issues related to it are one of the larger error sources influ- 
encing global modeling and its comparison to direct mea- 
surements. Note that many data products derived from di- 
rect measurements use a charge exchange cross section as 
well, as part of the ionization channels acting on neutrals. 
Therefore, the discussion of this systematic error applies 
to these deriv e d dat a pro ducts as well. For exam ples see 
iBzowski et all (|2008D and iRichardson et all (|2008( ) in this 
issue, where the uncertainty about the charge exchange 
cross sections is echoed in the interpretation of the pickup 
ion results, or the solar wind slowdown results, respectively, 
in terms of the interstellar H density. 



followed as a separate plasma fluid which interacts with the 
main solar wind protons. 

This list is not comprehensive, but is meant as a sample 
of the type of issues that are outstanding for the business 
of global heliospheric modeling, nonwithstanding past and 
present progress on multiple fronts (numerous citations are 
omitted here for the sake of brevity). Further progress, as 
well as extensions and refinements of additional lines of 
model physics, will improve the realism of all the models 
over time. 



5. Conclusions 

We investigate in this paper global heliospheric 
plasma/neutral models from five groups, first the plasma 
parts by themselves, then the fully self-consistent models. 
For the latter, the neutral species are modeled with either 
one, three, or four fluids, or on a particle-kinetic level. 
Performing model runs with exactly the same boundary 
parameters and physics included, we arrive at the following 
conclusions. 



4.3. Additional physics 

Finally, there are systematic errors influencing global helio- 
spheric modeling whose magnitude is difficult to assess or 
sometimes not yet explored. Neglecting interplanetary and 
interstellar magnetic fields, for example, excludes a whole 
suite of possible heliospheric asymmetries, shifts in the he- 
liospheric boundaries, and influences on the neutral hydro- 
gen distribution even in the inner heliosphere. The pres- 
ence of magnetic fields also typically allows for temperature 
anisotropics and turbulence in the plasma, and there is ev- 
idence that the solar wind (plasma) velocity distribution 
is non-Maxwellian already at a 1 AU distance, which most 
global models do not yet address. Further away from the 
Sun, the proton distribution functions are driven away from 
equilibrium by the effects of charge exchange, which calls 
for a fully kinetic plasma - neutral gas numerical modeling 
strategy eventually. 

The 3D, time-dependent solar wind in real-time differs 
from what most models currently feed into their simula- 
tions. Similarly, the solar irradiance depends non-trivially 
on time and on position in the heliosphere. Additional sim- 
plifying assumptions often made include the restriction of 
the particle species to electrons, protons, and neutral hy- 
drogen, and omitting heavier ion species (including alpha 
particles) in the solar wind, and heavier particles in the 
interstellar medium. The influence of high-energy particles 
such as cosmic rays (anomalous and galactic) should be 
taken into account, however their effect on the heliospheric 
structure is most likely not significant. 

The pickup process, i.e., the dynamical process of ac- 
celerating a newly born ion into the plasma bulk flow, also 
is often not handled in sufficient detail in global models. 
Many times, global modeling assumes instantaneous pickup 
for simplicity. It would be more realistic to fairly treat the 
pickup ion evolution, accounting for plasma-wave or tur- 
bulent interaction, and in general accounting for the non- 
Maxwellian character of the pickup ion distribution. A first 
level of refinement is taken in the Sch model (Fa hr et al.l 
I2OOOI) as used in this paper (section 3), where pickup ions 
are not absorbed instantaneously into the main plasma, but 



1. Although very different numerical strategies and ap- 
proximations have been chosen for the five heliosphere 
models presented in this paper, the results all qualita- 
tively agree. In many respects, even the models' quanti- 
tative predictions are similar. They agree in particular 
about the location of the upwind termination shock and 
the upwind heliopause. The discrepancies for termina- 
tion shock and heliopause in the five investigated models 
range from a few percent in the nose direction (<7%) to 
<14% in the downwind direction. The upwind distance 
of the bow shock disagrees by up to 15%. Also largely 
independent of the modeling strategy in this sense is 
the velocity of neutrals entering region 4 through the 
upwind termination shock (<11%). 

2. The pileup of neutral H in the hydrogen wall is sensi- 
tive to the modeling strategy, and the maximum density 
of neutral H differs by about 60% between the extreme 
cases. The column density through the hydrogen wall 
does not seem to vary; the hydrogen wall is either steep 
and narrow or small and broad. The strength of the in- 
terstellar bow shock and the associated post-shock ve- 
locity is the driver for the height of the hydrogen wall, 
and the same mechanism leads to a variation in the fll- 
tration, with the smaller hydrogen wall generally leading 
to less filtration (larger neutral density entering through 
the termination shock). The neutral H distribution in 
the inner heliosphere is therefore moderately sensitive 
to the strength of the interstellar bow shock. This is re- 
markable as the bow shock is at the farthest heliospheric 
distance. 

3. The strength of the bow shock also anticorrelates with 
its resulting distance from the sun. In comparing the 
five modeling strategies, the bow shock strengths dif- 
fer by 90% between the extremes of the five models. 
The bow shock is strongest for a two-fluid model, and 
turns out progressively weaker if neutrals are modeled 
with more and more fluids, and is weakest in the kinetic 
models. This behavior influences the neutral results in a 
systematic way, with the filtration being the strongest 
in the four-fiuid case, weaker for the five-fluid model, 
and weakest in the kinetic models. There are exceptions, 
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however, as the two-fluid model presented here yields a 
filtration as weak as the kinetic models. 

4. There are two discernible reasons for the different bow 
shock strengths in the models. First, different numerical 
strategies are used to model the bow shock itself, rang- 
ing from shock-capture methods to smeared-out discon- 
tinuities. Second, depending on the neutral modeling 
strategy, different amounts of secondary neutral hydro- 
gen make it to the region upstream of the bow shock, 
and these neutrals have a much larger kinetic tempera- 
ture when modeled with particle-kinetic methods rather 
than fluids. Both these factors prime the interstellar 
plasma through charge exchange upstream of the bow 
shock, and weaken the bow shock. 

5. Global heliospheric models without neutrals (section 2) 
do not reproduce many of the salient plasma density, 
velocity, and temperature features of the heliosphere ev- 
idenced by models with self-consistent neutrals (section 
3), as can be seen, for example, by comparing Figure 
[T^ to Figure [3^. Naturally, the absence of the contribu- 
tion of ISM neutral ram pressure to the pressure balance 
leads to an enlarged heliosphere in the plasma-only case 
(Table[2]vs. Tabled . This hence underlines the fact that 
the inclusion of neutrals in a global heliosphere model 
- in any self-consistent way - is critical to achieving 
physically meaningful results. 

The fluid and kinetic neutral atom models agree to 
within about 10% in some quantifiable measures, such as 
the location of the principal heliospheric discontinuities, 
which is similar to the uncertainties due to numerical al- 
gorithms, cross sections, grid size and resolution. Larger 
differences of about 50% exist in the details of hydrogen 
distribution function (hydrogen wall magnitude and neutral 
velocity distribution functions) between the models based 
on kinetic and hydrodynamic neutral atom descriptions. 
The uncertainties in our knowledge of the interstellar con- 
ditions, of charge exchange cross sections, inclusion of the 
MHD effects missing from the models discussed, and the un- 
explored effects of additional physics are further expected 
to modify the results by a similar amount. 
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